%%  Streamlines/stub
% Budget file
B=readBud([basename,'.bgt']);

% Compute stream function values
B=mf_Psi(B);                dPsi=0.2;   Prange=ContourRange(B,dPsi,[],'Psi');

% streamlines(o,ax,Psi,prange,varargin)

load underneath  % to get gr object

% set figure
fig=figure();                       % create handle
Ax1=subplot(1,1,1);
plottop;

for i = 1:length(B)
    if  i == 1
        [sliceStreamH]    = streamslice(gr.XM,gr.YM,gr.ZM,B(1).Psi,[gr.xm(1) gr.xm(end)],[gr.ym(1) gr.ym(end)],[],Prange);
        set(sliceH,'FaceColor','interp','FaceLighting','gouraud','EdgeColor','none','DiffuseStrength',2)
    else
        cdata = sliceCData(gr.XM,gr.YM,gr.ZM,C(i).values,[gr.xm(1) gr.xm(end)],[gr.ym(1) gr.ym(end)],[],hrange);
        for ii = 1:length(sliceH)            
            set(sliceH(ii),'CData',cdata{ii})
        end
    end
    drawnow;
end

%%  plot cocentration contours
C=readMT3D('MT3D001.UCN');  hrange=ContourRange(C,10); % get contour intervals
load underneath  % to get gr object

% set figure
fig=figure();                       % create handle
Ax1=subplot(1,1,1);
plottop;

% slicecontour loop
for i = 1:length(C)
    if  i == 1
        [sliceH]    = slice(gr.XM,gr.YM,gr.ZM,C(1).values,[gr.xm(1) gr.xm(end)],[gr.ym(1) gr.ym(end)],[],hrange);
        set(sliceH,'FaceColor','interp','FaceLighting','gouraud','EdgeColor','none','DiffuseStrength',2)
    else
        cdata = sliceCData(gr.XM,gr.YM,gr.ZM,C(i).values,[gr.xm(1) gr.xm(end)],[gr.ym(1) gr.ym(end)],[],hrange);
        for ii = 1:length(sliceH)
            set(sliceH(ii),'CData',cdata{ii})
        end
    end
    drawnow;
end
